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ABSTRACT 

We show that dark matter substructure in galaxy-scale halos perturbs the time delays between images in 
strong gravitational lens systems. The variance of the effect depends on the subhalo mass function, scaling as 
the product of the substructure mass fraction and a characteristic mass of subhalos (namely (m^) / (m)). Time 
delay perturbations therefore complement gravitational lens flux ratio anomalies and astrometric perturbations 
by measuring a different moment of the subhalo mass function. Unlike flux ratio anomalies, "time delay 
millilensing" is unaffected by dust extinction or stellar microlensing in the lens galaxy. Furthermore, we show 
that time delay ratios are immune to the radial profile degeneracy that usually plagues lens modeling. We lay 
out a mathematical theory of time delay perturbations and find it to be tractable and attractive. We predict that 
in "cusp" lenses with close triplets of images, substructure may change the arrival-time order of the images 
(compared with smooth models). We discuss the possibiUty that this effect has already been observed in RX 
J1131-123L 

Subject headings: dark matter — gravitational lensing — time 



1. INTRODUCTION 

The cold dark matter (CDM) paradigm predicts that 5%- 
10% of each galaxy's mass is bound up in subhalos left 
over from the hierarchical formation process. In the Local 
Group, the predicted number of dark matter subhalos signif- 
icantly ex ceeds the o bserved number of dwarf galax y satel- 
lites (e . g..lMooreetal. 1999; Klypin et al. 1999; Strigari et all 
l2007al: iKoposov et al.l |2008|) . The abundance of substruc- 
ture depends on competition between accretion of new subha- 
los from the environ ment and destruction of old subhalos by 
tidal forces (e.g., Tavlor& Babul 2001, 2004; Benson et ajj 
2002; Zentner & Bullock 2003; Kou shiappas et al 2005 
Oguri & Leejl2004t Ivan den Bosch et al. 2005; Z entner et all 



Also, the number of subhalos that "light up" and be- 
come visible as satellite galaxies depends on whether subha- 
los are able to retain their gas against photoevaporation, and 
on the efficiency of galaxy formation in low-ma ss systems 
(e.g. [Buflock et al. 2000; Somerville 2002; Kr avtsov et all 
l2004t IKoposov et aL 2009: Maccio et al. 2009) . Measuring 
the amount of substructure in galaxy halos, and how it varies 
with galaxy mass, environment, and redshift, therefore pro- 
vides unique access to the astrophysics of galaxy formation 
on small scales. 

We still know very little about the physical properties of 
the dark matter particle, but a number of specific models 
have been proposed: dark matter could be sterile neutrinos, 
or supersymmetric particles, or a manifestation of extra di- 
mensions, or even a product from the decay of any of these 
particles (e . g..lDodelson & Widrowl lT994t ICheng et al.ll2002l: 
iF'engI 120051: IStripari et al.ll2007bl) . AU of those possibilities 
are compatible with observations that probe the universe on 
large scales. However, they make some different predictions 
about the amount of dark matter substructure; for example, 
any type of "warm" dark matt er can lead to a suppression 
of po \ yer on small scales (e.g.. ICohn etal.ll2000t iDave et all 
120011; IZentner & BuUockll2003b . Studying galaxy substruc- 
ture provides the opportunity to test such models and obtain 
important astrophysical evidence about the fundamental na- 



ture of dark matter 

Strong gravitational lensing is a simple geometric phe- 
nomenon that gives us a valuable tool for studying mass in 
distant galaxies, including substructure. Lensing effects are 
succinctly encoded in the properties of the time delay surface. 
Consider light emitted by a background source at angular po- 
sition u from the center of a foreground lens with (projected) 
potential 0. ' If the light reaches us from the angular image po- 
sition X, the excess travel time relative to a hypothetical light 
ray that travels directly from the source with no deflection is 



t{x) = to 



1 2 

-\x-u\ -(j)(x) 



c Dis 



By Fermat's principle, images form at stationary points of 
the time delay surface, i.e., positions x such that Vt{x) = 0. 
(This directly yields the lens equation, usually written as 
u = x-W4>.) The magnifications of the images are deter- 
mined by second derivatives of r (and hence depend on sec- 
ond derivatives of the lens potential cfi). What is usually 
measured is the differential time delay between two images, 
Atji = Tjxj) - T (xi). (See the review of strong lensing by 
iKochanek et 1111200 4 for more details.) 

Strong lensing provides the only way to detect sub- 
structure directly (i.e., by virtue of its gravity) in galax- 
ies outside the Local Group. Small mass clumps in the 
lens galaxy can strongly perturb lensed images. The spa- 
tial perturbations are determined by first derivatives of 
the lens potential: they have angular scales of milli- 
arcseconds for " millilensing" by dark matter subhalos (e.g., 
Mao & Schneider 1998; Metcalf & Madau 2001; Chiba 200^ 



Dalai & Kochaneb i2002i: iMetcalf et al.i i2004c iChen et all 
2007 ), or micro-arcsecond s for " microlensing" b y stars 
in the lens galaxy (e.g., |W itt et al. 1995; Wvithe et alj 
120001: ISchechter & Wambsganssii2002i; iKochanek et al.ii200'7r 

' The lens potential </> is a scaled version of the two-dimensional gravita- 
tional potential. Specifically, it satisfies the Poisson equation V'</< = 2E/Ecrit 
where S is the surface mass density of the lens, and Scrii is the critical surface 
density for lensing. 
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TABLE 1 
Reference IMAGES: "Fold" 



TABLE 2 
Reference IMAGES: "Cusp" 
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Note. — Reference smooth images for our sample fold lens 
similai' to PG 1115+080. All coordinates are given in arcseconds 
with respect to the center of the lens galaxy at (0,0). The time 
delays are given with respect to the leading image. The sign of 
the magnification (p) represents the image parity. Col. 1 gives the 
traditional image names for PG 1115+080, while Col. 2 gives our 
new labels: "M" and "S" indicate (respectively) an image at a local 
minimum or saddlepoint of the time delay surface, while "1" and 
"2" indicate (respectively) the leading and trailing image of each 
type. 



ICongdon et alj l2007l) . Since lensing magnifications de- 
pend on second derivatives of the lens potential, they are 
even more sensitive to substructure: magnification pertur- 
bations can be of order unity and ar e therefore very appar- 
ent, especially in 4-imag e lens es (e.g ., iMetcalf & Zhaoll2002l: 
iKeeton. Gaudi & Pettersil200ll2005h . 

At optical and X-ray wavelengths quasar emission regions 
are small enough that lens flux ratios are sensitive to both 
dark matter s ubhalo s and stars (e.g., jK ochanek et al. 2007; 
Keeton et al] l2006t iBlackburne et all .2006.: Poolev et al.. 



20071) : this makes it difficult to isolate milhlensing and study 
dark matter substructure. By contrast, at radio wavelengths 
the quasar source is thought to be large enough to smooth over 
the effects of stars and b e insensitive to microlensing (but see 
iKoopmans & de Bruvnll2000.) . The flux ratios of radio lenses 
(especially 4-unage, or quad, lenses) have therefore been the 
tool of choice for studying millilensing. The amount of sub- 
structure needed to explain radio lens flux ratios is broadl y 
consistent with CDM predictions ( Dalai & KochaneklllO^ . 
The greatest limitation of this met hod is the small num ber of 
radio quads known currently (e.g.. lBrowne et al1l2003h . 

In this paper we show that lens time delays provide an ex- 
citing new way to probe dark matter substructure, with sev- 
eral distinct advantages. As we shall see, time delays are not 
affected by microlensing, so we can use optical as well as 
radio data to probe dark matter substructure; this is impor- 
tant because there will be thousands of lenses discov ered in 
new o ptical surveys (e.g., Fassnacht et al. 2004; Kuhl en et alJ 
12004: Ma rshall et al) 120051: iKochanek et alj I2OO60. comple- 
menti ng new samples of radio lenses (e.g., iKoopmans et alj 
l2004l) . Time delays are sensitive to the mass function of dark 
matter subhalos in a way that flux ratios are not; so they offer 
a good opportunity to probe the masses of subhalos in distant 
galaxies. By all indications the theory of time delay millilens- 
ing is very tractable; having a formal theory will provide a rig- 
orous foundation for substructure studies, and may even allow 
us to do some of the statistics analytically. Appropriate time 
delay measurements are feasible now, and truly revolutionary 
datasets will become available in the foreseeable future. 

For the purpose of concreteness, we assume a cosmology 
with Qm = 0.3, = 0.7, and Hq = 70 km s"' Mpc"', and we 
adopt specific values for the lens and source redshifts (moti- 
vated by particular lens systems, as discussed below). Mod- 
ifying the cosmology or the lens and source redshifts would 
simply rescale the time delays through the lensing time scale 
to- 



Note. — Similar to Table[T] but for our sample cusp lens similar 
to RX Jl 131-1231. We again introduce new image labels such that 
Ml is the leading minimum, M2 is the trailing minimum, SI is the 
leading saddle, and S2 is the trailing saddle. 

2. METHODS 

2.L Reference smooth lens 

We examine the effects of substructure by comparing the 
properties of images produced by a galaxy with substructure 
to those produced by an equivalent smooth galaxy. We model 
the smooth galaxy using a pseudo-Jaffe profile with scaled 
surface mass density 
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where Ecrit = {c^Ds)/{4ttGDiDis) is the critical surface den- 
sity for lensing, where Di, D^, and D;, are angular diame- 
ter distances between the observer and lens, the observer and 
source, and the lens and source, respectively. We have writ- 
ten the equation for a circular lens, but it is straightforward 
to obtain an elliptical model by replacing r with an appropri- 
ate elliptical radius. The pseudo-Jaffe model is equivalent to 
the singular isothermal model for r <C a, but the density falls 
more quickly at large radii r S> a to keep the total mass finite, 
Mtot = TTfl/jEcrit. The scale radius a can therefore be thought 
of as a sort of truncation radius, although the truncation is not 
sharp. We set a to be 300 kpc, although the particular choice 
has little effect on our results provided it is well outside the 
Einstein radius (see Section 12.41 for more discussion). The 
parameter /^tot is the Einstein radius of the lens. 

We calibrate the pseudo-Jaffe model by fitting it to observed 
lens systems. The fit is not quantitatively precise because 
we model the lens environment using a simple external tidal 
shear, when in truth the lens may lie in a group of galaxies 
that has a more comph cated effect (e.g.. lKundic et alj|1997l : 
iMomcheva et al.ll2006l) . For the present work, we seek only 
to obtain reasonable smooth mass models, not to replicate any 
observed lens in rich detail. 

We create a lens with a "fold" image con figuration simi- 
lar to that of PG 11 15H-080 (IWeymann et aUPigSO) . The lens 
galaxy redshift is z/ = 0.3 1 and the source redshift is Zs = 1 72, 
so the lensing time scale is to = 46.5 days (see Equation ([T]i). 
Our reference smooth model is circular with Einstein radius 
bioi = l."16, and it has external shear 7 = 0.12 at position an- 
gle 0^ = 65° (east of north). We place a source at position 
(-0."032,0."1 18) relative to the center of the galaxy, and ob- 
tain the four images listed in Table [T] and shown in Figure [T] 
below. 

We also create a lens with a "cusp" image configuration 
similar to that of RX J1131-1231 dSluse et al.. ,2003.) . The 
lens redshift is zi = 0.295 and the source redshift is Zj = 0.658, 
yielding a lensing time scale fo = 45.5 days. Our refer- 
ence smooth model has Einstein radius /jtot = 1"85, ellipticity 
e = 0.16 at position angle 6*^ = -56° (east of north), and ex- 
ternal shear 7 = 0.12 at position angle 6^ = -83°. We place a 
source at position (-0."549,-0."142) relative to the center of 
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Fig. 1 . — Upper panels: the squares show the images for our sample fold configuration (compare with Table[T). In terms of the traditional image labeling 
for PG 1115+080, we have M1=C, M2=Ai, S1=A2, and S2=B. Circles show the Einstein radii of mass clumps; panel (a) includes a key with m = 10** Mq for 
reference. The three panels show random realizations of a substructure population in which the substructure mass fraction is /, = 0.01, and the subhalo mass 
function is dN/dm oc over the range IO'-IO'Mq. Lower panels: sample fight curves for image M2. In these examples, we assume image Ml brightens 
by 20% at f = and returns to its original flux at f = 1 day. We plot the magnitude difference Am = mmi ~"1mi ■ In each panel the dotted curves shows the fight 
curve for the reference smooth model, while the solid curve shows the actual light curve for the specific realization of substructure. Notice that the substructure 
perturbs both the flux ratio and time delay of image M2 relative to image Ml. (The results for other images are similar.) 



the galaxy, and obtain the four images Hsted in Table |2] and 
shown in Figure |9] While there are quantitative differences 
in the time delays between the fold and cusp lens examples, 
many of the scalings and conceptual conclusions remain the 
same. Therefore in most of the presentation we focus on the 
fold lens for clarity. The exceptions are discussed in Section 



We should comment on our labeling of the lensed images. 
Images are traditionally labeled by letters (e.g.. A, B, C, D) 
using criteria chosen by the discoverers; different criteria have 
been used, with the result th at existing labels are not very in- 
formative. [S^afc^^Uiami (|2003) instead advocate labeling 
the images 1-2-3^ in arrival-time order, arguing that this la- 
beling is robust (at least for smooth mass distributions; see 
Section [331 l. and presenting rules for using the image config- 
uration to determine the labels. We assert that it is useful to 
have labels convey not only the time ordering but also the type 
of image as well. We therefore use "M" or "S" to indicate (re- 
spectively) an image at a local minimum or saddlepoint of the 
time delay surface, combined with "1" or "2" to indicate (re- 
spectively) the leading or trailing image of each type. Thus, 
the four images arranged in arrival-time order are Ml, M2, 
SI, and S2. In Tables [T] and |2] we give the correspondence 
between our labeling scheme and the traditional letter labels 
forPG 1115H-080andRXJ1131-1231. 

2.2. Subhalo population 

While detailed semi-analytic substructure models a re now 
available (e.g., Tavlor & Babul 2001, 2004; B ensori et al.l 
20021: Zentner & Bullockl 120031; Koushiappas et7l.l l2004t 



the average surface mass density in subhalos at radius r is 
i^s{>') = fs i^iotir), where fs is the substructure mass fraction. 

CDM simulations predict that the subhalo mass function 
is ap proximately a power law, dN/dm oc with /3 sa -1.8 
(e.g.. iGhigna et alJl2000l: iHelmi et aUlMl iGao et alJIIOOl 
iDiemand et alJ2007h . In order to explore how subhalo masses 
affect time delays, we consider masses in some finite range 
nil <m<m2. The mean subhalo mass is 



1 + /3 



i+li i+f) 
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'2 "'1 

Equivalently, we may write the lower mass limit in terms of 
the mean mass and dynamic range q = ni2/ni\ as 

, , 2 + /3 q^^l^-l 



nil 



(4) 



' 1+/3 ^2+/3_l ■ 

Another quantity that will prove to be useful is the mean 
squared mass. 



(m^) = (m) 
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Qguri & Lee ,2004; van den Bosch et al., .2005; .Zentner et alj 
2005h . for pedagogical purposes it is attractive to adopt sim- 



1+/3 3+/3 V^'"" ' ' ' 
When q = I the mass function reduces to a Dirac (5-function 
and all subhalos have the same mass. 

Figure [T] shows three realizations of a subhalo population 
in which the substructure mass fraction is fs = 0.01 and the 
subhalos have masses between nii = 10^ and m2 = IO^Mq. 
For reference, this mass function has (m) = 6.2 x IO^Mq and 
{m^)/{m)=2.8 x IO^Mq. 

2.3. Subhalo models 
For a subhalo of mass m, it is useful to define a scaled mass 



pie but reasonable assumptions about the mass function and 
spatial distribution of subhalos (e.g., Metcalf & Madau 200 H 
IChibal 120021: iDalal & KochaiieM [20021: iMetcalf et alJ l2004).2 
We assume that subhalos trace the total mass distribution, so 

^ A few millilensing studies have worked directly with A'-body simu- 
lations, but both the original simul ations and the len s ing calculations are 
computationally ex pensive tBradac et alj I2002L I2004t lAmara et ai] 120061 : 
IMaccio et aU2006l) . 
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which has dimensions of area. Here /?Ein is the Einstein radius 
of the subhalo if it is a point mass. In most of our analysis 
we do treat subhalos as point masses for simplicity. The lens 
potential of a subhalo is then 



) = — In r , 



(7) 
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while the deflection angle is 



d(j) 
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We also consider subhalos modeled as truncated isothermal 
spheres. The truncation occurs in the three-dimensional den- 
sity profile, so we have p oc out to the truncation radius r,, 
and p = for r>r,. The lensing deflection angle works out to 
be 



r < r, : a = ■ 



r > r, : a= — 
nr 



1-1 



1/2 



(9) 



It is straightforward to obtain the potential by integrating, (p = 
J adr, but it is not instructive to write the (somewhat long) 
expression here. 

2.4. Calculations 

To obtain a mass model with substructure, we leave a frac- 
tion (l-fs) of the dark matter mass in a smooth component, 
and replace a fraction fs with subhalos. The expected number 
of subhalos is 



f.Mtot Trabtotfs 



(10) 



and we draw the actual number from a Poisson distribution 
with this mean. The subhalo positions are assigned randomly 
as follows. For the pseudo-Jaffe model the cumulative proba- 
bility distribution for the radius is 



Pr(r)-- 
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We pick a random value for Pr uniformly between and 1 , 
and then invert Equation (fTTT i to find the random radius. The 
explicit inversion is 



Pr(2-Pr) 
2(1 -Pr) 



(12) 



(Again, these equations are for a circular lens, but the exten- 
sion to an elliptical mass distribution is straightforward.) We 
then pick a random azimuthal angle uniformly between and 
In. 

The pseudo-Jaffe model formally extends to infinity, but we 
do not actually need to consider subhalos at large radii be- 
cause they have little effect on the time delays. In the Ap- 
pendix, Equation (IA13I ) allows us to calculate the rms error 
that we make on time delays if we neglect subhalos beyond 
some radius Rq. We use this equation to set the threshold 
radius to ensure that our time delay errors are less than 0.001 
days. The fact that time delays are not very sensitive to distant 
subhalos is the reason that our results do not depend strongly 
on the choice of the pseudo-Jaffe "truncation" radius. 

The subhalo masses are assigned randomly from the power 
law mass function. The cumulative probability distribution 
for the mass is 



Pmim) = 



i+g 1+/3 



1+/3 



(13) 



We pick a random value for f ,„ uniformly between and 1 , 
and invert Equation dTsT l to find the random mass. (When the 



mass function is a (5-function, all halos have the same mass 
and this process is unnecessary.) 

The number of subhalos in a simulation varies between a 
few hundred and a few hundred thousand, depending on the 
subhalo masses and the substructure mass fraction. We use a 
tree algorithm ( Barnes & Hu t 1986) to compute the net lens- 
ing potential, deflection, and magnification quickly and effi- 
ciently. All of this analysis is included in an updated version 
of the public software package gmvlens (Keeton 2001). 

To find the images as they are affected by substructure, 
we note that the shifts in the image positions are small (e.g., 
IChen et ani2007h . so it is efficient to start from the smooth 
images and perturb the positions iteratively: 
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where a*'' and /i*'^ are the deflection vector and magnifica- 
tion tensor, respectively, evaluated at the current position .t*'\ 
(Here / labels the iteration step, not the image index.) Note 
that a;*'' - a*'' represents the source position associated with 
the current image position, so (5u*'' is the offset in the source 
plane. When the image position correctly solves the lens 
equation, vanishes and the iteration process converges. 

A massive clump near one of the "macro-images" can in 
principle split it into several "milli-images".^ We check for 
this possibility by comparing the parity of the recovered im- 
age against the parity of the reference image. If we start with a 
minimum but recover a saddle (or vice versa), we know there 
are extra images and we have found one. Note that if we start 
with a minimum and recover a minimum (or likewise for a 
saddle), that does not prove there are no extra images. How- 
ever, we use the minimum/saddle cases to estimate the frac- 
tion of simulations that yield extra images, and find it to be 
small. We throw away the small number of simulations that 
yield extra images. 

Once we have found the image positions, we use Equa- 
tion ([T]i to compute the differential time delay between image 
pairs, Af/j = T(xj)-T{xi). 

3. RESULTS 
3.1. Time delay perturbations 

We begin by examining our sample fold lens in the pres- 
ence of substructure, as shown in Figure[T] In the bottom row 
of that figure we show sample light curves for image M2 (as- 
suming for pedagogical purposes a simple square wave varia- 
tion in the source); the results for the other images are similar. 
The substructure clearly changes both the flux ratio and time 
delay between the images. In these examples, the flux ratio 
is perturbed by 0.1-0.2 m ag and the tim e delay is perturbed 
by some fraction of a day. lOguril (l2007h also showed (in less 
detail) that substructure can affect lens time delays, but he fo- 
cused on this as a source of noise in lensing measurements of 
the Hubble constant. We instead propose that time delay per- 
turbations provide a new way to detect and study substructure 
itself, if lens time delays can be measured well enough. 

The remainder of the paper is devoted to exploring how 
time delay perturbations depend on the subhalo population. 
To that end, we start with a simple case in which all the sub- 
halos are point masses with m = IO^Mq, and the substruc- 

' With point mass clumps, or indeed any clump with a cuspy central den- 
sity steeper than S oc r"' , there is formally a "micro-image" very close to 
each clump. However, these images are highly demagnified and hence unim- 
portant. They are absent when the clump central density is shallower than 
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Fig. 2. — Histograms of the time delays between tlie images, for 10* 
Monte Carlo simulations witli substructure mass fraction /j = 0.01 and sub- 
halo mass m = 10^ Mq. (The vertical axis scale is arbitrary.) The vertical 
lines indicate the time delays for the reference smooth model. Superposed on 
each histogram is a Gaussian (dotted line) with the same mean and variance 
as the simulated data. 
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Fig. 3. — Histograms of the time delays between different images (rows), 
for different values of the subhalo mass (columns). We now plot 5{At), the 
difference between the time delay with substructure and the delay for the 
reference smooth model. All cases have lO** Monte Carlo simulations with a 
substructure mass fraction /, = 0.01. Notice that the horizontal scale varies 
from one column to the next. (The vertical axis scale varies as well, but 
the units are arbitrary anyway.) Equivalent Gaussians are again shown with 
dotted hnes. 

ture mass fraction is = 0.01. For each substructure real- 
ization we compute the time delays between the images; we 
then repeat the process 10"* times and make histograms of 
the time delays, as shown in Figure |2] For each image pair, 
the average time delay matches the prediction of the refer- 
ence smooth model very well, which makes sense because the 
smooth model can be thought of as an average over the sub- 
halo population. These three time delay histograms appear to 
be Gaussian. (The histogram of the time delay between the 
close images M2 and S 1 is somewhat different, as discussed 
in Section [331 ) We attribute the Gaussianity to the Central 
Limit Theorem: the lens potential is a sum of many random 
terms, and such a sum tends toward a Gaussian distribution if 
the mean and variance of each term are finite, which we ver- 
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scaling ct; oc (/^m)'/^ from Equation lAlOt in the Appendix (after we fit for 
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ify in the Appendix.^ The standard deviation is ct, = 0.24-0.26 
days for all three image pairs shown in Figure |2] 

3.2. Subhalo mass and substructure mass fraction 

Next we consider the effects of changing the subhalo mass 
(Figure O and the substructure mass fraction (Figure |4|l. In 
order to place different image pairs on a common scale, we 
plot histograms of the time delay perturbation S(At), defined 
to be the difference between the time delay with substructure 
and the time delay for the reference smooth model. All of the 
histograms appear to be nicely Gaussian, to be centered on 
zero, and to have widths that increase with both the subhalo 
mass and the substructure mass fraction. 

We can quantify this last point by plotting the time delay 
scatter a, as a function of the subhalo mass, for different val- 
ues of the substructure mass fraction, as shown in Figure|5] In 
the Appendix we use the Central Limit Theorem to make an 
initial prediction of the scaling of the time delay scatter (see 
Equation ( lAlOb ): 

a,«(/,m)'/2. (15) 

To test this prediction, we use the results from simulations 
with /; = 0.01 to fit for the proportionality constant, and then 
plot the analytic curves alongside the simulation results in 
Figure |5] The agreement is striking given the simplicity of 

* The Central Limit Theorem has not been used in millilensing before, 
because it does not obviously apply to flux ratios. See the Appendix for more 
discussion. 
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the analytic prediction. Strictly speaking, the analytic scaling 
appears to overestimate the time delay scatter when the sub- 
halo mass is large. We can identify at least two possible ex- 
planations. First, when /, is fixed increasing the subhalo mass 
decreases the number of subhalos, which may make the Cen- 
tral Limit Theorem less applicable and invalidate an approx- 
imation in Equation iA6i . Second, as the subhalo mass in- 
creases the perturbations to the image positions also increase, 
which may invalidate one of our starting assumptions in the 
Appendix. We are working on a more sophisticated theory 
of time delay perturbations to address both issues. For now, 
though, we consider the success of the initial analytic predic- 
tion to be very encouraging. 

The scaling of the time delay scatter with subhalo mass has 
an important corollary: there is no measurable effect on time 
delays from low-mass objects, such as stars. This means we 
can study millilensing at any wavelength without "contami- 
nation" from microlensing. 

3.3. Subhalo mass function 

So far we have assumed that all subhalos have the same 
mass; now we consider a power law mass function. Figure |6] 
shows the time delay scatter as a function of the dynamic 
range q = m2/mi and power law slope (3 of the mass func- 
tion. The dynamic range has a significant effect: the broader 
the mass function, the more scatter there is in lens time de- 
lays. The power law slope affects the time delay scatter as 
well, but much more modestly. For comparison we also plot 
the analytic scaling from Equation ( IA9I ) in the Appendix, 



1/2 



(16) 



where the ratio (w^) / (m) can be found from Equation (|5]l. 
The analytic formula does not match the simulation results 
perfectly, but it seems remarkably good for something so sim- 
ple. We take this as an indication that the theory of time delay 
millilensing will be analytically tractable. Having a rigorous 
theory of millilensing is possible only with time delays; an- 
alytic results for flux ratio millilensing are elusive because 
flux ratios are highly nonlinear and the Central Limit Theo- 
rem does not apply (see the Appendix). 

3.4. Internal structure of subhalos 

So far we have treated the mass clumps as point masses, 
partly for simplicity and partly because we know the gravity 
outside any spherical mass clump is the same as that of a point 
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mass. However, dark matter subhalos in general have some 
spatial extent, and if they overlap lensed images this may 
be important for millilensing. In this pilot study we mainly 
want to determine when the point mass approximation is rea- 
sonable and when we need to worry about subhalo size and 
structure. We therefore consider a simple model of isother- 
mal spheres that are truncated (in three dimensions) at some 
radius r,. Figure [T] shows how the time delay scatter changes 
when we vary the truncation radius while keeping the subhalo 
mass fixed. (The key scale is the dimensionless truncation ra- 
dius, r, = rf/REin, where /?Ein is the Einstein radius of a point 
mass with the same mass, as defined by Equation 

Depending on the mass, subhalos need to extend beyond 10 
or more Einstein radii before the internal structure becomes 
important. The obvious next step is to model tidal truncation 
and estimate the sizes of realistic subhalos. There are some 
important effects to consider: tidal forces vary with position in 
the parent halo; and tidal truncation also corresponds to mass 
loss (so, strictly speaking, it is not right to vary the truncation 
radius while keeping the subhalo mass fixed, although that is 
a useful pedagogical exercise). These effects are incorporated 
into semi-analytic substructure models, so in follow-up work 
we will study time delay millilensing using realistic subhalo 
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Fig. 8. — Histogram.s of the time delay of the close image pair M2/S1 for 
our sample fold lens (compare Figs.[2]and[3). The substructure mass fraction 
is /, = 0.01, and the subhalo mass is indicated in each panel. The vertical 
hnes indicate the time delay for the reference smooth model. Superposed on 
each histogram is a Gaussian (dotted line) with the same mean and variance 
as the simulated data. 

populations drawn from those models. 

3.5. Arrival-time ordering of tensed images 

ISaha & Williams! (I2003h suggested that the arrival-time or- 
dering of lensed images is a robust prediction of smooth lens 
models (at least, those with relatively simple angular struc- 
ture). With the help of a few simple rules, they noted, it is 
usually easy to identify the two minima and two saddlepoints 
in a quad lens, and then to deduce the arrival-time ordering. If 
the situation is not immediately obvious, a simple lens model 
makes it clear. We have therefore chosen to use the image 
classification and ordering as the basis of our labeling scheme: 
in arrival-time order we have the leading minimum (Ml), the 
trailing minimum (M2), then the leading saddle (SI), and fi- 
nally the trailing saddle (S2). 

In the fold image configuration we have considered, sub- 
structure does not affect the arrival-time ordering. The rea- 
son is that the time delay perturbations are small compared 
with most of the time delays themselves. It might seem that 
the close image pair M2/S 1 could provide an exception, be- 
cause the smooth model time delay (0.16 day) is shorter than 
some of the substructure perturbations we have seen. How- 
ever, by definition a saddlepoint in the time delay surface must 
be higher than a nearby local minimum, so image S 1 must trail 
M2 and the overall ordering must be preserved. One conse- 
quence is that the histogram for the M2/S1 time delay can- 
not remain Gaussian as the substructure perturbations grow, 
as shown in Figure [8] This implies that the substructure in- 
fluences both images at some level, although the /S.tM2.s\ his- 
tograms are broad enough that it seems fair to say that sub- 
structure affects even nearby images quite differently. The 
only way this general analysis can be violated is if a massive 
and concentrated clump lies close enough to images M2 and 
SI to change the global topology of the time delay surface; 
but such a situation would be immediately apparent from the 
image configuration (see Saha & William&„2003) . 

The issue of arrival-time ordering is strikingly different for 



a cusp image configuration like that in RX J1131-1231. Fig- 
ure 13 a) shows the images along with selected arrival-time 
contours for our smooth mass model. In this example, and 
indeed in all s mooth models of RX J1131-1231 examined by 
[Morgan et al.l ((2006), the image labeled Ml is the leading im- 
age. However, as shown in Figure |9f b), substructure can re- 
verse the arrival-time order of the images Ml and M2, such 
that the overall ordering is M2/M1/S1/S2. In terms of the 
topology of the time delay surface, we can still say the sad- 
dlepoint must be higher than the two minima on either side, 
so SI must trail both Ml and M2. However, there are no 
such restrictions on the heights (or depths, equivalently) of 
the two minima with respect to each other. It is fairly easy for 
substructure to modify one or both of the minima enough to 
change their relative heights, and hence the arrival-time order- 
ing of images Ml and M2. In fact, some 27% of realizations 
with modest substructure (substructure mass fraction /, = 0.01 
and subhalo mass m = IQ^Mq) predict an arrival-time reversal 
(see Figure [TOb. 

We believe image reversal by substructure may already 
have been observe d in RX J1131-1231 , based on the time 
delays reported by [Morgan et al.l (l2006l) . The time order- 
ing observed among the three bright images is M2/M1/S1^ 
with AfM2,Mi = 2.20!j;^| days, AfM2,si = ll-98!|:p days, and 
AfMi.si = 9.61^j 57 days. These time delays present two puz- 
zles. First, image M2 is observed to lead image Ml, whereas 
smooth models predict the opposite. Second, the observed 
M2/S1 and Ml/Sl time delays are an ord er of magnitude 
longer than predicted by smooth models. Mor gan et al.l (l2006h 
focused their attention on the second problem, and showed 
that placing a single, massive (> 5 x 10"^ M©) subhalo near 
image SI could explain the long M2/S1 time delay. Their 
models still predicted Ml to be the leading image, though. 

As we have demonstrated, a population of subhalos might 
be able to reverse the ordering to make M2 the leading image. 
It will still take a lot of modeling to draw quantitative con- 
clusions about the substructure required to explain the image 
ordering and the long time delays; we will present the mod- 
eling details elsewhere (C. R. Keeton & L. A. Moustakas, in 
preparation). For now, we mainly want to introduce the idea 
that the temporal ordering of images in a cusp configuration 
can be changed by substructure — and that this effect may al- 
ready have been observed. 

3.6. Macromodel uncertainties 

Having shown that substructure affects lens time delays, 
we still need to consider how well the perturbations could be 
detected in observed lenses. The question is whether uncer- 
tainties in the "macromodel" used to characterize the smooth 
component of the lens mass distribution could mask the ef- 
fects of substructure. One key concern is the radial profile 
degeneracy: especially in quad lenses, varying the radial den- 
sity profile of the lens galaxy can rescale the time delays while 
leaving o ther observables unchanged (e.g., Falco et al. 19851 
lKeeton&Kochaneklll997l: ISahiilOOOt iKoch anek 2002). At 
first glance, this would seem to make changes induced by sub- 
structure degenerate with changes in the radial profile. 

However, the radial profile rescaling affects all images in 
the same way, while substructure affects each image differ- 
ently. So with more than one time delay (as in a quad lens), 
ratios of the delays cancel the rescaling and isolate substruc- 

^ We have switched to our labeling scheme using the identifications M1=C, 
M2=B, S1=A, andS2=D. 
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Fig. 9. — (a) Image configuration for our sample cusp lens, assuming a smooth mass distribution (compare with Table|2j- In terms of the traditional image 
labeling for RX J1131-1231, we have M1=C, M2=B, S1=A, and S2=D. Positions are measured with respect to the center of the lens galaxy at (0,0). Selected 
contours of constant arrival time are also shown. Counting contours shows that the minimum containing image M2 is shallower than the minimum containing 
image Ml, hence M2 trails Ml. (b) Example of the same lens with substructure (with fs = 0.01 and m = 10** Mq). Now it is apparent that the minimum containing 
the image labeled M2 is actually deeper than the minimum containing Ml . In other words, substructure has reversed the airival-time ordering of these two images 
such that M2 actually precedes Ml. The ordering of SI and S2 remains unchanged. (Note that the contour levels are different in the two panels; they are chosen 
to best reveal the minima and saddlepoint in each case.) 
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Fig. 10. — Histograms of the time delays between the close images in our 
example cusp lens. The images labeled Ml and M2 lie at minima of the time 
delay surface, while SI lies at a saddlepoint. In smooth models the image 
ordering is always M1/M2/S1/S2. However, in Monte Carlo simulations with 
a substructure mass fraction fs - 0.01 and subhalo mass m = 10^ Mq, about 
27% of the realizations have AfMi,M2 < which means image M2 actually 
leads image Ml. In such cases the overall ordering in M2/M1/S1/S2. As 
usual, vertical lines indicate the time delays for the reference smooth model, 
and a comparison Gaussian is superposed on each histogram. 

ture effects, as shown in Figure[TT] The radial profile degener- 
acy will therefore not be a fundamental limitation in attempts 
to probe substructure with time delays. This very promising 
inference is currently empirical, but we expect it can be made 
rigorous with a full theory of time delay millilensing. Since 
radial profile effects do not cancel perfectly (there are small 
differences in the widths of the time delay ratio histograms in 
FigurefTTTi. the theory will also be useful in understanding and 
correcting for these residual effects. 

Another source of uncertainty is the angular structure 
of the macromodel. While it is common to fit ob- 



served lenses with an ellipsoidal mass distribution and ex- 
ternal shear (and we have used such a model to con- 
struct our mock lenses), real galaxies may be more com- 
plex. In particular, observed (e.g.. Bender e t al.l Il989t : 
ISaglia, Bender & Dressier" 1993'; Rest et al. 2001) and simu- 
lated (e.g., Hevl, Hernauist & Scerael 1994; Naab & Burker] 
120031; IJesseitetal.1 120051: ,Naab et all 120061) elliptical galax- 
ies often contain mild departures from elUptical symmetry 
that can be described with ot = 3 and m = 4 multipole terms 
in the surface mass density. lEvans & Witt! (l2()03i) suggested 
that such multipole terms might provide an alternative to 
clumpy substructure when fitting observed lenses. Although 
multipoles turn out not to provide a satisfactory e xplanation 
of observed flux ratio anomalies ( Kochanek & Dalall 120041; 
ICongdon & Keetonll2005l: lYoo et al.ll2005i. j2006i) . it is impor- 
tant to consider how they might influence the analysis of time 
delays. 

We can do this by taking our mock lenses generated with 
substructure and trying to fit them with multipole models. For 
this exercise we adopt 0."003 errorbars for the image posi- 
tions, which is typical of observations with the Hubble Space 
Telescope (HST) and radio interferometry, and 5% time de- 
lay uncertainties, which is similar to the preci sion achieve d 
for quad lenses today (see the compilation by IOgurill2007h . 
In general, we find that multipole models can formally pro- 
vide a good fit to the image positions and time delays, but 
many of the models have unreasonably large multipole terms 
that lead to unrealistic galaxy shapes, as shown in Figure [T2l 
To avoid the unreasonable models we can impose some pri- 
ors on the multipole terms. Since the multipole moments in 
observed and simulated elliptical galaxies are typically at the 
percent level, a simple first step is to adopt Gaussian priors 
on the multipole amplitudes and 04 with a mean of zero 
and standard deviation of 0.02. (In the present analysis we 
do not impose any constraints on the orientations of the mul- 
tipole terms.) With these mild amplitude priors, multipole 
models fail (i.e., they are inconsistent with the data at more 
than 95% confidence) for 27 of the 100 mock lenses we exam- 
ine. We are in the process (J. van Saders et al., in preparation) 
of deriving more sophsticated priors on multipole terms from 
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Fig. 11. — Examining the radial profile degeneracy in our mock fold lens. The smooth mass component has a power-law density profile p oc r or 
equivalently a surface density profile ft oc while the substructure is still given by our pseudo-Jaffe model. Varying the profile rescales the time delays (left 

and middle columns). However, the time delay ratio (right column) is largely unaffected: the mean ratio stays the same, while there is a small change in the 
scatter. (The time delays between images Ml, M2, and S2 are used for illustration; time delays involving image SI could be used as well.) Time delay ratios 
therefore allow us to probe dark matter substructure without worrying about the radial profile degeneracy. 




Fig. 12. — Lensing critical curves for multipole models fit to six random mock lenses. We fit the image positions and time delays using multipole terms that 
are unconstrained (dotted curves) or have Gaussian priors (solid curves; see text). We also consider adding the image flux ratios as constraints (dashed curves). 
Among these lenses, cases b and c are the most glaring examples of failed multipole models; cases d and f are less dramatic but failures nonetheless. 



a large sample of galaxies in the Gal axy Evolution fr om IVIor- 
phologies and SEDs survey (GEMS: lRix et al.l|200 4'). and we 
expect that those will be even more valuable in constraining 
macromodels and thus aiding the identification of time delay 
anomalies. 

It is perhaps artificial to insist that we must identify anoma- 
lous lenses on the basis of time delays alone, when flux ratios 
are also available to provide additional constraints. If we fit 
the flux ratios along with the image positions and time delays, 
we find that multipole models fail (at 95% confidence) for 52 
of 100 mock lenses even if we use generous 10% flux error- 
bars and do not impose any priors on the multipole terms. In- 
terestingly, only 28 of those failures would be identified on the 
basis of flux ratios alone; in 24 of the mock lenses, it is only 
when we fit the time delays and flux ratios simultaneously that 
the anomalies become clear We infer that time delays and 
flux ratios contain different and complementary information 
about the small-scale structure of the lens potential. 

While we plan to address the issues of multipole priors and 
complementarity of time delays and flux ratios in detail in 
future work, we offer the preliminary conclusion that macro- 



model uncertainties can be controlled well enough that we can 
expect to use time delays to learn about substructure. 

4. DISCUSSION AND CONCLUSIONS 

We have shown that substructure in lens galaxies modifies 
the time delays in multiply-imaged gravitational lenses. The 
amplitude of time delay perturbations depends on the abun- 
dance of substructure, the mass function of subhalos, and to 
some extent the internal structure of subhalos as well. This 
phenomenon, which we call time delay millilensing, has some 
very attractive properties. First, the effects of substructure on 
time delays are fairly easy to understand, both conceptually 
and quantitatively. Second, time delay perturbations are unaf- 
fected by stellar microlensing or extinction in the lens galaxy. 
Third, time delay ratios are immune to the radial profile de- 
generacy in lens modeling. 

Furthermore, time delay millilensing is sensitive to the 
mass function of subhalos in a different way than flux ra- 
tio and astrometric millilensing. For point mass perturbers, 
the lensing optical depth depends on j R\^^^(dN / dm)dm cx 
J m(dN /dm)dm. The potential of a clump scales with 
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its mass, so time delay perturbations depend in full on 
J np' (dN/dm)dm. (See the Appendix for more details.) 
For comparison, position perturbations have a scale of 
/?Ein C!C m'/^ while magnification perturbations are dimen- 
sionless, which means that astrometric millilensing de- 
pends on / ifp^^ {dN/dm)dm while flux ratio millilensing de- 
pendson f m(dN/dm)dm (compare Dalai & Kochanek 2002; 
Dobler & Keetonl 120061: IChen et alJ 120071: IShin & Evans! 
2008h . 

We therefore suggest that lens time delays should join im- 
age positions and flux ratios as tools for studying dark matter 
substructure in distant galaxies. Since the three observables 
measure different moments of the subhalo mass function, the 
combination of all three will ultimately provide a powerful 
way to test CDM predictions and even to explore the nature of 
dark matter itself, given that different dark matter candidates 
may produce different mass functions at subgalactic scales, as 
a function of lens galaxy mass and redshift. 

Before we try to use time delays to measure substructure, 
it is vital that we understand whether systematic uncertain- 
ties in lens models could corrupt substructure constraints. We 
have explicitly shown that the radial profile degeneracy is not 
a significant problem, as time delay ratios essentially fac- 
tor out the profile dependence. We have drawn the prelim- 
inary conclusion that uncertainties in the angular structure 
of the lens potential can be dealt with, using additional flux 
ratios data and/or independent constraints on the shapes of 
galaxy mass distributions. Other potential concerns may in- 
clude the environment of the lens galaxy (Keeton & Zabludo3 
2004. and (sub)s t rucfiire along the fi ne of sight (|Keetoni20()3l: 
ChenetalJ[200l iMetcal j l2005allbl: [ui^ |2008|). For all of 
these issues, the important point is again that non-local fea- 
tures in the lens potential affect the images in some coor- 
dinated way, whereas substructure affects the images differ- 
ently. 

We also need to predict time delay perturbations for 
subhalo populations that are more realistic than the sim- 
ple ones we have used (for pedagogical purposes) here. 
This work is underway, using existing semi-analytic sub- 
structure models. By incorporating important physi- 
cal effects such as accretion of new subhalos and tidal 
truncation and disruption of old subhalos, these models 
yield subhalo populations whose spatial distributions and 
mass functions agree well w ith results from A'-body sim- 
ulatio n s (e.g., | Tavlor & Babul 2001, 2004; Benson et al' 
20021; IZenti^& Bullock .2003: Koushiappas et al.. ,2004: 
Qguri & Led 120041: Ivan den Bosch eTall 120051: IZentner et alj 
20051) . Our follow-up calculations will offer greater insight 
into the subtleties of time delay perturbations in a broader 
range of realistic circumstances, and better guidance as to how 
well we need to measure time delays if we want to probe dark 



matter substructure. 

We have noted that in extraordinary lenses with close 
triplets of images in a cusp configuration, it may be possi- 
ble to detect the effects of substructure through unambigu- 
ous changes in the arrival-time ordering of the images (com- 
pared to expectations for standard macromodels). We ex- 
pect, though, that the real power of time delay millilens- 
ing will lie in its application to large ensembles of lenses 
with measured time delays. At present, time delay uncer- 
tainties are typically at the level of 1 day (see IOgurill2007l 
for a compilation). Improvements of a factor of a few are 
desired: our preliminary estimates indicate that improving 
the uncertainties to ^0.5% would yield constraints on sub- 
structure at the leve l of 0.3-0.5 dex in a given 4-image lens. 
ICoUev et alj ( l2003h have shown that coordinated global op- 
tical monitoring can yield time delays with uncertainties of 
less than 0.1 day. Monitoring at X-ray wavelengths, where 
variability can be rap id and dramat i c, can yield even more 
precise time delays: IChartas et alJ (l2004l) measure a time 
delay of 3.58 ±0.14 hours between the two close images 
in PG 1115H-080 with Chandra and XMM-Newton observa- 
tions. With some effort, such techniques can probably be 
extended to a modest sample of time delays in the near fu- 
ture. Over the longer term, the advent of large time-domain 
surveys is expected to yield thousands of lensed quasars 
with time delavs (e.g.. iFa ssnacht et al. 2004; Kuhlenetal] 
120041: [Marshall et alJl2005l: iKochanek et al] [2006). Even bet- 
ter would be a dedicated space platform capable of pre- 
cise monitoring of len ses over extended periods of time (see 
iMoustakas et al.l2008 ). Looking to the future, measuring pre- 
cise image positions, flux ratios, and time delays in a sample 
of ~100 quad lenses could yield substructure measurements 
at the 20%-30% level or better Since the measurements will 
span the redshift range 0.2 < z < 1, this will provide a unique 
opportunity to study not just the amount of dark matter sub- 
structure but its evolution as well. 
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APPENDIX 

STATISTICS OF TIME DELAY PERTURBATIONS 

In this Appendix we present an initial theory of the statistics of time delay perturbations. This theory rests on three assumptions: 
(1) subhalos are statistically independent (of each other; the probability of having a subhalo may still depend on position); (2) 
subhalos may be treated as point masses; and (3) time delay perturbations are dominated by changes in the lens potential (as 
opposed to changes in the image positions). The first assumption means that all subhalo positions and masses are drawn from 
the same probability distribution. This simplifies the theory, and seems reasonable for lensing since we work in projection and 
subhalos that are near each other on the sky are most likely well separated in space. The second and third assumptions will be 
relaxed in future work, but are adequate for our present goal of understanding in a general way how the time delay scatter depends 
on the amount of substructure and the subhalo mass function. 

The lens time delay has the form given by Equation By assumption #3, if we want to understand how substructure affects 
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Afi2 = t{xi)-t{x2), it is sufficient to study the potential difference A(j)i2 = (f){xi)-(f>{x2)- Any change in the time delay is simply 
given by the change in the potential difference, multiplied by the time scale to. 

The potential due to a collection of subhalos is the sum of contributions from individual subhalos. If the subhalos are point 
masses, the potential difference has the form 

N , I I 

A012 = > ip(ri,mi), where ^p(r,m) = —In- -. (Al) 

^ TT |2:2-r| 

If we can show that the mean and variance of (p are finite, then we can use the Central Limit Theorem to infer that the sum A0i2 
will be approximately Gaussian, and to determine the variance of A(j)i2. 

In order to compute the mean and variance of (p, we need to average over a subhalo's position and mass. Let Prir) be the 
probability distribution for the subhalo's position, while pmim) = (1 /N)dN /dm is the probability distribution for its mass. Then 
the average of any function / has the form 

(/) = / dm p,„(m) [ dr pAr) f . (A2) 



It is useful to note that smoothed substrucure density distribution, K.s(r), may be written in terms of the position probability 
distribution: 



: y <i»J ^ vr/?|i„ pr(r) = J dm ^ m pAr) = N (m) pr(r) . (A3) 



Now in order to determine the mean and variance of ip, we need to compute the moments 

(i^) = / dm pm(m) I dr pAr) — In ■ 



TT \X2-r\ 

1 f la;i-rl 
— / dr nAr) In f-! \ , (A4) 

m^ ( \x\ — r\ 



Nt: J ' \x2 — r\ 



(</5^) = j dmp^(m) j dr pAr) ^ ^In 



\x2-r\ 

J_ (4') ...... r,J-i-H^' 

Ntt^ (m 



dr KAr) I InT^—- 7 I • (A5) 



In both cases we used Equation (IA3b to substitute for pAi") in terms of k^. We also identified the mass integrals as yielding the 
averages (m) and {m^). The variance is then 

At the second step, we note that {(fi^) ^ 1 /N while (ip)^ ^ 1 /N^, so when there are many subhalos we may neglect the second 
term. 

The integrands in Equations (IA4I) and ( IA5I ) diverge near xi and X2, but only logarithmically, and such a divergence is still 
integrable. As a result, as long as Kj(r) does not diverge badly near the image positions, the integrals are well behaved, so the 
mean and variance of tp are finite. This means the Central Limit Theorem holds, and we can infer that the potential difference 
A012, and hence the time delay Afn, will have a distribution that is approximately Gaussian. (This argument also applies to 
clumps that are not point masses, because then any divergence in p will be softer than logarithmic.) 

Furthermore, thanks to the Central Limit Theorem we may compute the variance in A(j)i2 as the quadrature sum of the variances 
of all the individual ip terms — or, in this case, a simple multiplication by A^: 

al^Nal^\^^fdr.Ar)(ln\{^y. (A7) 
^ TT^ [m) J \ \x2-r\J 

And of course a, = focr^. Heuristically, we may summarize this result as 

a,aU.,^j , (A8) 

where is indicative of the amount of substructure. While this expression oversimplifies the dependence on the spatial dis- 
tribution of subhalos, it is conceptually instructive. In our simulations, we write nAr) = fs Ktoi(r), and we keep the total mass 
distribution Ktot fixed while we change the substructure mass fraction Therefore it is formally correct to write 

/.W ■ (A9) 
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This is our general conclusion for how the time delay scatter depends on the amount of substructure and the subhalo mass 
function. Note that if all subhalos have the same mass m, this simplifies to 

a,<x(f,my/\ (AlO) 

It is instructive to consider whether the Central Limit Theorem can be applied to flux ratios as well as time delays. For point 
mass clumps, the magnification of an image at position x has the form 

[(l-^..)(l-<f>yy)-4>l]~\ (All) 

where 

m (x-Xif-(y-yif 



:=-E 



TT \x-Xi\'^ 



(A12) 



and there are similar expressions for 0, , and (/)^, . Now when we try to compute the mean and variance of each term b y integ rating 
over Xi, the integrand diverges non-integrably near x, and so the mean and variance diverge (also see Witt et al. 1995). The 
Central Limit Theorem therefore does not apply. The magnification probability distribution p(fj,) may still be perfectly well 
defined, but it need not be Gaussian. Inde ed, p(m) can be computed analyt ically for a uniform distribution of point masses, and 
it is distinctly non-Gaussian (see §11.2 of lSchneider. Ehlers & Falcolll992 ). There may be circumstances (say, certain types of 
spatially extended clumps) for which the Central Limit Theorem does apply to magnifications, but the theorem is not universally 
applicable as it is for time delays. 

Returning to time delays, there is one last bit of theory that is useful. We intuitively expect that clumps far from the images have 
little effect on the time delay, but we can quantify this statement. Consider for simplicity images at xi = (d,Q) and X2 = (-d,Q), 
and let us examine only subhalos at radii r > Rq^ d. We return to Equation iAli . plug in for xi and X2, and then take a Taylor 
series expansion to lowest order in d/r. 



4 « — ^ / de drr K^r) Ad^ « \^ dr . (A13) 



TT^ (m 



Notice that we have changed notation and written to highlight that this is an error term — specifically, the rms error we make by 
neglecting subhalos beyond Rq. We have computed the error in the potential, but we can relate it to the rms error in the time delay 
via Cf = foe^. (We have verified Equation (IA13b using direct simulations of subhalo populations extending to large radii.) We 
have assumed for simplicity that the substructure mass distribution is circularly symmetric, n^ir) = Ks(r), but the generalization 
to elliptical symmetry is straightforward. For any reasonable substructure distribution, Ks(r) decreases as r — > cx), so the integral 
converges. Not only that, but e, is a monotonically decreasing function of R{). In other words, the farther the subhalos are from 
the images, the less effect they have on the time delays. If we can tolerate some small error in time delays, we may neglect 
all subhalos beyond some threshold radius ^o- Once we specify the time delay error tolerance, we can solve Equation ( |A13l l to 
compute the threshold radius Rq. 
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